Deformation of inherent structures to detect long-range correlations in supercooled 
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We propose deformations of inherent structures as a suitable tool for detecting structural changes 
underlying the onset of cooperativity in supercooled liquids. The non-affine displacement (NAD) 
field resulting from the applied deformation shows characteristic differences between the high tem- 
perature liquid and supercooled state, that are typically observed in dynamic quantities. The 
average magnitude of the NAD is very sensitive to temperature changes in the supercooled regime 
and is found to be strongly correlated with the inherent structure energy. In addition, the NAD 
field is characterized by a correlation length that increases upon lowering the temperature towards 
the supercooled regime. 

PACS numbers: 61.43.Fs,64.70.Q-,05.20.Jj 
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I. INTRODUCTION 

Many liquids, when quenched fast enough, enter a su- 
percooled regime that is signaled by an enormous increase 
in viscosity eventually leading to glass transition. Several 
experimental as well as theoretical studies have shown an 
accompanying qualitative change of the dynamics: While 
individual particle dynamics dominate at high tempera- 
tures, particle motions become increasingly more coop- 
erative and heterogeneous as the temperature is lowered 
EHIJ]. From the sizes of cooperatively rearranging re- 
gions, a growing dynamical length scale can be extracted 
near the glass transition 

One of the most puzzling features of the supercooled 
regime and the glass transition is the apparent lack 
of structural changes underlying the dramatic slowing 
down of the dynamics. Some progress has recently been 
achieved in linking structural and dynamical changes 
near the glass transition. For a two-dimensional model 
system, a connection between dynamical heterogeneities 
and local crystalline order was suggested in Such 
a connection appears to be absent in a different two- 
dimensional system, where instead dynamical hetero- 
geneities seem to be correlated with local, short-time dy- 
namics 0] and localized soft modes @ . A somewhat dif- 
ferent path has been followed by several researches, try- 
ing to relate dynamical properties to inherent structures, 
which are the local minima of the collective potential 
energy It has been shown that inherent struc- 

tur es g overn mechanical properties of amorphous solids 
[14l4l6| . For supercooled liquids, more and more evidence 
has been collected that qualitative changes in the dynam- 
ical behavior are accompanied by changes in the inherent 
structures 

EES El- The local diffusivity, for example, 
was suggested to be related to the basin depth of local 
inherent structures [l9j ]. Similarly, unstable modes of 
saddle points in the energy surface were found to be of 
major importance, not only for the mobilities of particle 
clusters E] but for the glass transition in general [2l| . 



Collective rearrangements in the inherent structures cor- 
responding to the dynamics of supercooled liquids are 
reported in El- Specially designed simulations probing 
point-to-set correlations that are inspired by random first 
order theory have revealed evidence for growing amor- 
phous order in the supercooled regime (23j . In some very 
recent works, external forcing is applied to (some parts 
of) the system in order to extract characteristic sizes of 
cooperative regions using mode-coupling El or density 
functional theory El- Under shear flow, it was found 
numerically that coop erative, mobile regions can form 
anisotropic bands |2q . 

We have found striking similarities in the onset of co- 
operative behavior in dynamics and in the non-affine part 
of the inherent structure response to external deforma- 
tions]^ . That approach was motivated by a recent the- 
ory [28j based on a general framework of non-equilibrium 
thermodynamics. This thermodynamic treatment re- 
quires information about the system's response to an ap- 
plied, static deformation and suggests that the reversible 
part of glassy dynamics changes considerably when ap- 
proaching the glass transition. Above the glass transi- 
tion, the particles can follow an imposed deformation 
more or less freely, whereas closer to the glass transition 
the particle movement becomes a hopping-like transition 
between different basins of attraction of the underlying 
inherent structures Ej El ■ By imposing static deforma- 
tions of inherent structure configurations, we observed in- 
deed a profound difference in the NADs when approach- 
ing the glass transition. From a systematic finite-size 
scaling analysis, we found that the NAD field is charac- 
terized by a static correlation length, that is growing as 
in usual critical phenomena Ej E] . This length detects 
growing structural correlations underlying the growing 
dynamical length scale obtained from particle dynamics 

aim- 

Here we present an extensive study of the NAD field 
introduced in Ej Ei that further demonstrates how 
this quantity can effectively detect the onset of the co- 
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operative dynamics and shows a great potential to probe 
the underlying structural correlations. The distribution 
of lengths changes from exponential to power law upon 
entering the supercooled regime, and we can rational- 
ize such change in terms of a crossover from a viscous 
liquid to a regime dominated by elastic effects. The 
mean displacement length depends exponentially on the 
inherent structure energy, as we also discuss using a sim- 
ple toy model, and confirms the existence of two well- 
distinguished regimes as a function of the temperature. 
We analyze different measures of correlations in the di- 
rection of the NAD field and discuss their analogies with 
observations in the dynamics . Finally we use a coarse- 
graining procedure to extract the characteristic size of 
correlated regions observed in the snapshots: we discuss 
different definition of this length-scale and we perform a 
finite-size scaling analysis over different model systems, 
confirming the critical-like behavior at low temperatures. 

The paper is organized as follows. The model glass for- 
mers used and the numerical simulations are described in 
Sect. In] In Sect. IIII1 we briefly review the method pro- 
posed in 27] to extract displacements of inherent struc- 
tures. Our numerical results for the NAD are presented 
in Sects. HVl and IVl Correlations between these displace- 
ments are analyzed in Sect. |V] We focus in Sec. [IV] on 
characterizing the NAD lengths and we introduce a sim- 
ple model to rationalize the results. In Section IVll we ex- 
tensively describe the coarse graining procedure for the 
analysis of the NADs introduced in Ref. [30[ and ap- 
ply it to the different model systems to extract the tem- 
perature and system size dependence of the correlation 
length. Sect. |VlT] contains further discussion and conclu- 
sions. 



II. MODEL DESCRIPTION 
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50:50 


« 0.11 
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TABLE I. Some parameter values for the two models used. 
Tk denotes the (extrapolated) Kauzmann temperature and 
Tmct is the mode-coupling temperature. The mode coupling 
temperature is estimated in Ref. [34| and [35| for the BMLJ 
and BMSS models, respectively. For the BMLJ model, T K 
has been estimated numerically as Tk ~ 0.30 [3tj There 
have been several theoretical and numerical estimates for Tk 
in the BMSS model; 0.11 < T K < 0.14 have been reported 
in Ref. [33, [3S|. Here we take the numerical estimate given 
in Ref. [38|]. The onset temperature T a of the non-Arrhenius 
behavior of transport properties is estimated to be T a fts 1.0 
for the BMLJ model. T a depends somewhat on the quantities 
investigated and is to be taken as a rough estimate. A careful 
investigation of T a for the BMLJ model is presented in [iol ]. 
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do is the unit of length and density is chosen to 
be p — N/V = Uq 3 ■ A smooth cut-off is used by set- 
ting the potential to $ a b(r) = B ab (a - r cut ) 3 + C ab for 
a > r > ?* cu t = v3 and 3> a &(f) = Cab for r > a. The val- 
ues of B a b, C a b and a are fixed by imposing the continuity 
up to second derivative for $ a t(r). 



B a b = 



169e Ogb a 2 

r 3 > J ' 
' cut ' cut 



5e (Tab si2 
Cab = j , 

13 r cu t 



a = (15/13)r cut . 
The potential is then shifted to ensure that $ a b(a) = 0. 



The studies presented in this paper are based on com- 
puter simulations of two different, well-established mod- 
els for structural glasses (see e.g. HU). 

The first model is the three-dimensional, binary mix- 
ture of Lennard- Jones (BMLJ) particles introduced by 
Kob and Andersen [34[. Both particles have unit mass 
m and interact via a Lennard- Jones potential, $>ab(r) = 
^abii^ab/r) 12 — (cr aD /r) 6 ], where the diameters and in- 
teraction energies are given by uaa = &o, &ab = 0.8cto, 
(Jbb = 0.88(7o, £aa = £) ^ab = l-5e, and ebb = 0.5e. We 
have chosen cto as the unit of length and e as the unit 
of energy. The potential is truncated and shifted to en- 
sure $ aD (r C ut ) = 0, where the cut-off distance is chosen 
as r cut = 2.5a a b- 

The second model considered in this study is a 50:50 
binary mixture of soft spheres (BMSS) in three dimen- 
sions [23j . Both particle types have unit mass and the 
interaction potential is given by & a b{r) = e(cr aD /r) 12 , 
Cab = Ca + Ob, a,b = {A, B} , where the sizes of particles 
o a b are fixed by setting oa/ob = 1-2 and the effective di- 
ameter to one; that is (2a J 4) 3 + (2cr B ) 3 + 1(oa + o\b) 3 = 



Main parameter values for the two models are shown in 
Table |U For ease of notation, we use the same symbols 
for original and reduced quantities. 

For both model systems studied here, we have carefully 
prepared 10-100 independent samples for each temper- 
ature by molecular dynamics simulations starting from 
statistically independent, random initial configurations. 
Periodic boundary conditions were used in all cases. We 
have studied systems with number of particles N varying 
from 2000 up to 64000, but most of the data here refer 
to systems with N = 8000. The simulations were per- 
formed with the molecular dynamics simulation package 
LAMMPS ; 1 j j . The initial particle configurations were 
equilibrated at several decreasing values of temperature 
in the range 0.20 < T < 1.0 (BMSS), and 0.42 < T < 3.5 
(BMLJ). Comparing these temperature intervals with 
the mode-coupling temperature Tmct given in table U 
our simulations cover the high temperature regime down 
to the supercooled state and extend below the mode- 
coupling temperature. For all systems, slowly cooling the 
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configurations towards low temperatures was achieved by 
coupling the system to Nose- Hoover thermostat with pre- 
scribed, slowly decreasing temperature protocol. When 
the target temperature was reached, the temperature of 
the thermostat was held constant and the system was 
equilibrated at this temperature still in contact with the 
thermostat. We verified that no significant drift in the in- 
ternal energy nor in the kinetic temperature could be ob- 
served after the thermostat was switched off and the sys- 
tem was further evolved, now in the microcanonical en- 
semble. To test the equilibration of samples, we also com- 
pared the kinetic (Tkin) with the configurational temper- 
ature (T con f) (HI- The kinetic temperature is defined by 
NfksTkin = YljLi mv j i w he rc Nj denotes the number of 
degrees of freedom. While the kinetic temperature is en- 
tirely determined by the particle velocities, the configura- 
tional temperature T con f depends on the positions of the 
particles via fc B T conf = (| £\ Fi| 2 )/(- £\ V, • Fj). Here, 
Fj denotes the total force on particle j and V 3 = d/drj. 
We verified that kinetic and configurational temperature 
agree for our equilibrated samples within numerical un- 
certainties. Finally, we calculated two-time correlation 
functions and verified that no significant aging was ob- 
served in the equilibrated samples for a waiting time of 
the order of 40r Q , where r a is the structural relaxation 
time. 



III. METHOD FOR EXTRACTING NAD 

To extract non-affinc particle displacements we pro- 
ceeded as proposed in our earlier work [53]. In order to 
make the paper self-contained and settle the notation, we 
briefly review this method here. 

A. Afflne deformations 

We apply static, affine shear deformations to the par- 
ticle configuration by mapping the particle positions 
Yi — > rf , with rf = r.; +jyi& x , where 7 denotes the defor- 
mation amplitude. Since we aim at the non-affine part 
of the inherent structure deformations [28|] , we suggested 
in [27| the following procedure, summarized schemati- 
cally in Fig. [TJ start with configuration X = {r^} at a 
given temperature T . Prepare one configuration X dq by 
first applying the static deformation —> rf and sub- 
sequently finding the inherent structure corresponding 
to this deformed configuration. The other configuration 
X qd is prepared by first finding the inherent structure 
corresponding to the initial configuration X and after 
that subjecting the inherent structure configuration to 
the same deformation. Almost all subsequent analysis 
is based on the comparison between the configurations 
X dq and X qd , which we denote as non-affine displace- 
ment (NAD), dj = r dq — r qd . Thereby, we focus on the 
dependence of the NAD on the temperature T of the 
initial configuration and the amplitude 7 of the applied 



deformation. We ensure that the total displacement van- 
ishes, = 0, since a rigid translation can always be 
added and does not contribute to the NAD. 



B. Inherent structure generation 

From the equilibrated samples, we generate inherent 
structure (IS) configurations by minimizin g th e potential 
energy via the conjugate gradient method [43j. The min- 
imization is stopped when the potential energy change is 
less than a tolerance value 10 _7 e. We verified that the 
results are insensitive to a further decrease of the tol- 
erance level and that the mean inherent structure ener- 
gies for different temperatures agree well with published 
data [13, IH, [HI • For generating the inherent structure of 
deformed configurations, A dq , the minimization is per- 
formed in a deformed simulation box or, equivalently, 
using Lees-Edwards boundary conditions. 

C. Inherent structure properties 

Inherent structures can be characterized by their mean 
energy eis- In the inset of Fig. [5]we plot eis of the inher- 
ent structure X q as a function of T. Starting from a high 
temperature plateau, eis decreases upon cooling in the 
so-called landscape dominated regime [12]. In addition, 
we have performed various types of structural analysis 
on the IS of different model systems, using for example 
pair correlation functions, coordination numbers or bond 
order parameters, to characterize the differences between 
the particle configurations in the initial state and in the 
IS. As an example, in Fig. [2] we plot the averaged bond 
orientational order parameter (Qe) for the model BMLJ. 
(Qe) is calculated by averaging, for each temperature, 
over 100 independent samples of N = 8000 particles, the 
i-th order bond orientational order parameter [4oj ] with 
1 = 6. For particle j, this quantity is calculated as 

« = i/inh ( t WJ 2 ) ' m 

\m=— I / 

where Q\ m is the locally averaged bond orientational or- 
der parameter of order I and degree m as is defined in 
[H [13]. In the main frame of Fig. [2j (Q 6 ) is plotted 
as a function of temperature for the starting configura- 
tion X and its inherent structure X q . We observe that, 
in the range of temperatures where eis displays a strong 
dependence on T, the 6th order averaged bond orienta- 
tional order clearly increases for both X and X q . Also, 
as expected, the value obtained for X q is larger. In 
spite of this, the values of (Qe) are consistently much 
lower than those expected for various crystal lattices 
(e.g. Qe ~ 0.66 for perfect icosahedral ordering), demon- 
strating that there is no significant crystallinity present 
in the X and X q configurations [48j |. 
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IV. CHARACTERIZING THE NAD FIELD 

We have applied the procedure described in Sec. HTT1 
to well-equilibrated configurations of the two model sys- 
tems. Various temperatures of the initial configurations 
as well as different deformations are considered. Results 
are obtained as averages over independently generated 
configurations and error bars are calculated from the 
standard deviation. 



A. Distribution of NAD vectors 

We start by characterizing the distribution of the NAD 
vectors h(d). For the small deformations 10~ 5 < 7 < 
10~ 2 we are interested in, the distribution h(d) is found 
to be isotropic for homogeneous deformations and covers 
a range of displacement magnitudes. We therefore con- 
clude, that our regime of deformations is strong enough 
to allow the escape from local minima. On the other 
hand, the deformations are small enough such that no 
obvious trace of the imposed shear geometry is seen in 
h(d). This can be seen from Fig. [3J which shows the 
distribution of the Cartesian NAD components hi(\d a \), 
a = 1,2,3. The distributions h\ show a maximum at 
d a = 0. At high temperatures, /ii(|d Q |) decays expo- 
nentially, hi(\d a \) oc exp [— \d a \/d c ], with fit parameter 
d c w 0.1 of the order of the mean NAD length at this tem- 
perature. This behavior is reminiscent of the exponential 
distribution observed in Ref. [22j, also with 4 « 0.1, for 
displacements between inherent structures corresponding 
to the system's dynamics. 

This finding indicates again that the non-affine inher- 
ent structure deformations bear remarkable similarities 
to dynamical properties. At low temperatures, the dis- 
tribution hi is much narrower, in agreement with the 
impression from the snapshots shown in Fig. 2J More 
quantitatively, one finds that the distribution at low tem- 
peratures is no longer decaying exponentially but seems 
to exhibit power-law tails h\(\d a \) oc |d Q | _l/ with an ex- 
ponent approaching v w 2.5 at the lowest temperatures 
investigated, see the solid line in Fig. [3] It is interesting to 
note, that the power-law distribution /ii(|d a |) oc |do,|~ 5 / 2 
is predicted for the induced radial displacements in an 
elastic continuum around an expanded spherical shell 
[49l [50j . If we can identify the latter with localized re- 
arrangements, the qualitative change of the NAD distri- 
bution from exponential to a power-law at lower temper- 
atures apparently indicates the crossover from a viscous 
liquid to a regime with more pronounced elastic effects. 
From the goodness-of-fit, we determine the crossover 
temperature where the power-law and exponential distri- 
butions fit equally well to the data. We find the crossover 
temperature to be located in the interval T € (0.65, 0.70) 
for the BMLJ model (data not shown). 



B. Mean Mismatch length 

A first characterization of the NAD field {dj} is pro- 
vided by their mean length, 

d(T, 1 ) = {N-'J2tf) 1/2 : (2) 

3 

where the ensemble average is performed over the inde- 
pendently generated samples. Figure [5] shows d as a func- 
tion of temperature T for homogeneous deformations of 
different strengths 7 for the BMLJ model. Depending on 
the value of 7, we distinguish three different regimes. For 
small values of 7 (7 < 10~ 4 ), wc find that d is essentially 
independent of 7. This is an interesting observation as it 
suggests that we have reached a regime of weak pertur- 
bations, where details of the local deformation are less 
important for d. Moreover, d decreases drastically, by 
roughly two orders of magnitude, as the temperature is 
lowered from T > 1 to T ps 0.45. For larger deforma- 
tions 10~ 3 < 7 < 10 -1 , the decrease of J gets less and 
less pronounced, until d becomes essentially independent 
of temperature for deformations 7 > 0.1. The particu- 
lar form of the drastic decrease in d(T) depends on the 
model system. Qualitatively, however, the observations 
are the same also for the BMSS model. 

The decrease of d with temperature for small deforma- 
tions shown in Fig. [5] is quite similar to the decrease of 
the inherent structure energy eis (see [l3| and the inset 
of Fig. [S]). Therefore, we parametrically plot in Fig. [5] 
the average NAD length d versus eis for the different 
temperatures investigated. Our result indicate that d is 
indeed strongly coupled to the inherent structure energy. 
From Fig. [5] we observe that d increases exponentially 
with the inherent structure energy. More quantitatively, 
we observe a crossover from d ~ e _,3l ' eis ', /?i = 22.5 at 
low temperatures to d ~ e -fo\eis\^ p 2 _ g a ^ higher 
temperatures. Since the mean inherent structure ener- 
gies vary inversely proportional with temperature in the 
landscape- influenced regime [TH, we find that d ~ e~ A l T 
with a different constant A in the two regimes. The 
crossover from /3i to fa takes place at eis ~ —7.63, which 
corresponds to a temperature T ps 0.6, near to the inflec- 
tion point of eis(T), see the inset of Fig. [3J 



C. Simple toy model for inherent structure NAD 

Several features of the NAD field can be rationalized 
using a simple toy model. Let X denote again the state of 
our system and X q its closest minimum (inherent struc- 
ture). For simplicity, however, the toy model describes 
these states by a single scalar variable. The imposed de- 
formation of strength 7 deforms these configurations to 
X d = X + g and X qd = X q + g, respectively, where g 
denotes an average amount of displacement. The inher- 
ent structure of the deformed system is X dq = X q + nT, 
where T denotes a typical size of an inherent structure 
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basin. The integer n takes values n = 0,1,..., with n = 
corresponds to staying in the same minimum, n = 1 to 
the nearest minimum, etc. The NAD d is therefore given 
by d = nT — g and its mean squared average by 



d=[Y / Pn(nT-g) 2 } 



2ll/2 



(3) 



In Eq. ([3]), we have introduced the probabilities p n of 
finding X dq in the n-th nearest minimum. In the limit 
of Pn — > Sn,o, we find that d = <?, i.e. the NAD length is 
given by the imposed deformation. This is the case for 
extremely low temperatures and/or vanishing deforma- 
tions, not considered in the present study. For very large 
deformations compared to the basin size g > T, we ex- 
pect equal probability of entering a basin n in the neigh- 
borhood Am around m = g/T, where m 3> 1. Therefore, 
p n w 1/2 Am for \n — m\ < Am and zero else. Inserting 
this expression into Eq. ([3]) we find d « gAm/m AmT. 
Hence, in the case of very large deformations, the average 
NAD is expected to approach a limiting value. The simu- 
lation data in Fig. [5] seem to indicate such a trend. How- 
ever the mean length remains weakly sensitive to temper- 
ature variations at the largest deformation investigated. 

Finally, we consider the case of rather small deforma- 
tions, where one can approximately set p n = for n > 1. 
The important quantity is then p\, the probability of 
leaving the basin at the imposed deformation 7. Ac- 
counting in a rough way for the equilibrium distribution 
within the basin, we assume p\ = 7e~^ e , where e denotes 
a typical energy barrier and /3 = (fcBT) -1 . The values 
of 7 are restricted to 7 < 1 in order to ensure proper 
probabilities p\ > for all temperatures. With these 
assumptions, the mean NAD becomes 



d = r[ 7 (l - 2 7 )e 



-He 



■7 



211/2 



(4) 



where we have assumed g = , i.e. the important length 
for the imposed displacement is the typical basin size. 
Equation Q predicts a number of features that can be 
tested by simulations. First, for high temperatures, d 
reaches a limiting value d — > T^^(l — 7) independent of 
temperature, in agreement with simulation results (3l| . 
This limiting value decreases as 7 increases towards its 
maximum value. Second, for small deformations 7 < 1, 
Eq. (|4]) simplifies to d = T^/jexp [— /3e/2] and predicts 
a square-root increase of d with the strength 7 of im- 
posed deformation. In the simulations for the BMLJ 
model, d indeed increases with 7, the exponent, how- 
ever, changes with temperature from 0.2 at T = 1 to 1.0 
at T = 0.42 [3l|. Third, when the temperature is de- 
creased, Eq. ([5]) predicts a decrease of d according to the 
Boltzmann factor. This is in agreement with the results 
discussed in Sec. lIVBl (see also Fig|6|). Finally, Eq. <(4]) 
shows an exponential dependence of the average NAD 
on the barrier height e. If we can associate the inher- 
ent structure energy with an effective barrier height, this 
prediction is consistent with our numerical results (see 
Fig. More quantitatively, we assume that the effec- 
tive barriers are given by e(T) = — 2ft£;BT(eis(T) — eg), 



where eg is the high-temperature plateau of the inher- 
ent structure energy and the coefficient ft describes the 
dependence of d on eis in the low temperature regime, 
as introduced at the end of Section IIVBI A linear rela- 
tion between mean inherent structure energies and bar- 
rier heights was indeed found numerically for low temper- 
atures in the BMLJ model [5l| and similarly in a binary 
soft sphere model 2l|. For low temperatures, the mean 
inherent structure energy decreases as eig = e$ — o 2 /k^,T, 
where a 2 can be interpreted as the variance of the inher- 
ent structure ener gy d istribution in the Gaussian energy 
landscape model [l3j . In fact, we find in this regime 
that the height of effective barriers decrease linearly with 
temperature, e(T) — eo — ak^T for fceT < eo/a, where 
e = 2ft <r 2 and a = 2ft (e - eg). For the BMLJ model, 
we have eg w -7.55, e « -7.39, a 2 m 0.12, ft w 22.5, 
resulting in eo ~ 5.4 and a « 7.6. The value of the ex- 
trapolated barrier height eo is close to the one obtained 
from the mean waiting time for low- lying IS [5l| . 

Within the toy model Eq. @, the crossover from low 
to high temperature behavior observed in Fig. [5] could 
be interpreted as different regimes in the potential en- 
ergy landscape, with different relations between minima 
(inherent structures) and barrier heights. Interestingly, 
numerical simulations of randomly perturbed inherent 
structures for BMLJ model indicated that the number of 
saddles vanishes exponentially below T < 0.6 [52], close 
to the temperature corresponding to the inflection point 
shown in Fig. [5] as evaluated at the end of Section HVBI 



V. CORRELATIONS IN THE NAD FIELD 

Not only the length of the NAD vectors but also their 
spatial correlation change considerably between the liq- 
uid and supercooled regime. Spatial correlations in the 
NAD field can be clearly seen in the snapshots in Fig. [4] 
and were also reported in our earlier study [27| . In the 
following we thoroughly analyze them and use them to 
extract correlation lengths. 



A. Correlated mobility and directions 

We quantify correlations in the NAD field over a given 
distance r and for deformation amplitude 7 by 



C*(r, 7 ) 



m 2 ) 



(5a) 

C\\(r,7) = (52446(r-n j )), (5b) 
C±{r, 7 ) = (£di-df5(r-r ij )}, (5c) 



• ■j 



where r.^ denotes the distance between particles i and j 
and the average in the denominator is performed over all 
particles. The normalization is chosen such that Cg = 1, 
C11 = 1/3 and C± = 2/3 for perfectly correlated particles 
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and C a = for uncorrelated displacements, a = {5, IL-L}- 
The same correlation functions were used in Ref. Q in 
order to study correlated displacements in hard sphere 
colloids. The first function, Cs measures correlations in 
the mobility, i.e. the magnitude of the particle displace- 
ments (with the average subtracted), Sdj = \dj \ — (|d|), 
whereas C\^± is sensitive to correlated directions. The 
latter distinguish longitudinal and transverse correlations 
of directions, with df = d ; • rjj and d^ = dj — dfrjj, 
with f,-j the unit vector connecting particles i and j and 
dj = dj/|di|. 

Fixing the particle separation around the first neighbor 
shell and varying the deformation amplitude 7, we find 
that all three correlation functions ([5]) arc rather insensi- 
tive to 7 and start to decay only for 7 > 0.1 (not shown). 
This value is consistent with the treshold value deter- 



mined in Ref. [53| for shear-induced transitions between 
different inherent structures. Since we observe correlated 
rearrangements at smaller 7, wc would rather conclude 
that such large deformations induce transition between 
uncorrelated metabasins (see also Fig. IT2"j) . The behav- 
ior of Cy and C±, as well as the relation Cs > Cu + C 
found, are quite similar to the one observed in Ref. 
if our deformation 7 is identified with a properly scaled 
time (for not too short times). 

Next, we fix 7 to a typical small value (7 = 10~ 4 ), as 
already done in Sec. lIVI and study the spatial dependence 
of the correlation functions ([5]). We observe that the 
correlations in mobility measured by Cs increase as the 
temperature is lowered, but the spatial extent remains 
within a few particle diameters. Even shorter ranged are 
the transverse correlations C±. The longitudinal corre- 
lations Cu, however, increase significantly with decreas- 
ing temperatures as shown in Fig. [7J For small sepa- 
rations r, the correlation functions C a oscillate due to 
the short-range ordering as measured by the pair cor- 
relation function. For larger separations (r > 4), the 
oscillations have decayed and we observe an exponential 
decrease C\\ (r) oc e~ r ^^ . We find the correlation length 
£|| to increase from £|| sa 2.5 at high temperatures to 
almost 4 upon approaching the glass transition. Similar 
observations have been made in Ref. [54| in simulations of 
the NAD in amorphous solids and in Ref. @ for displace- 
ments in hard sphere colloidal systems when increasing 
the particle concentration. In the latter, however, the 
length scale was found to remain on the order of 3 par- 
ticle diameters, increasing significantly only beyond the 
glass transition. It was argued in Ref. [f| that the long- 
range correlations of the longitudinal displacements re- 
flect the string-like cooperative motion observed in com- 
puter simulations (5||. If this is indeed the case, our 
findings suggest that these string-like motions arise due 
to underlying string-like rearrangements between nearby 
inherent structures. At larger distances, one would ex- 
pect a hydrodynamic or elastic-like response where the 
correlations decay not exponentially but cx 1/r in three 
dimensions (56| . In agreement with Ref. @, we find no 
indication of such a behavior, probably because the corre- 



lation function has already decayed to such small values 
that the numerical data do not allow us to detect the 
expected power-law. 



VI. COARSE-GRAINED NAD 

We investigate correlations in the direction of the NAD 
vectors with the help of a coarse-grained displacement 
field around every particle, 



(6) 



Here, we have defined the orientations dj = dj/|dj| 
is the distance between particles j and k in 



and 



dq 



the inherent structure of deformed configuration. A dq , 



Nj = J2k Wb ( r 'jk) 1S th c number of neighbours of parti- 
cles j within a distance 6, and Wb(r) = 1 if r < b and 
zero else. Since the mean length of the NADs is strongly 
temperature-dependent (sec Fig. O, we use the normal- 
ized displacements in the definition ([5]) in order to sep- 
arate this aspect and focus instead on the correlations 
between vector orientations. 

When the coarse-graining length b is smaller than in- 
terparticle distances, only one particle contributes to the 
average in Eq. and Dj = dj. As b increases, more 
and more particles are involved in the average and the 
magnitude of Dj decreases. Figure [5] shows the function 
D(b) = (D 2 (6)) 1 / 2 obtained for the BMLJ model when 
subjected to homogeneous shear deformation with am- 
plitude 7 = 10~ 4 . The expected decrease of D with b is 
indeed observed. For a fixed distance b, D(b) increases 
monotonically with decreasing temperature. This behav- 
ior clearly indicates increasing correlations between par- 
ticle's NAD orientations as the temperature is lowered. 

Moreover, Fig. [5] shows that these growing correlations 
extend over distances larger than the particles diameter. 
Since D(b) measures the root-mean-square NAD direc- 
tion when averaged over a length b, we expect D to decay 
when b exceeds the size of a correlated region. Therefore, 
the decay D(b) exp [— b/£' D ] allows to define a length 
scale £' D over which the NAD vectors are correlated [lj| • 
At low temperatures, we find that D(b) is well described 
by an exponential decay over a wide range of smoothing 
lengths. The length scale £' D that we obtain from a least- 
square fit to D(b) is shown in the inset of Fig.[8]as a func- 
tion of temperature. At high temperatures (T > 0.8), 
an exponential decay of D (b) can be observed only in a 
narrow interval 2 < b < 4. In this regime, we find an 
approximately temperature-independent length scale of 
£' D ps 3. This length scale is somewhat larger than the 
size of correlated liquid structure measured by the first 
peaks in the pair correlation function. Since local rear- 
rangements have to include neighboring particles, it is 
not surprising that £' D is on the order of a few particle 
diameter in this regime. For decreasing temperatures, we 
observe an increase of the correlation length £' D to values 
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of about 6 particle diameters. Qualitatively, the same 
observations are made also for the BMSS model. 

In order to quantify the correlations in the NAD vec- 
tors in a different way and to show that considering the 
magnitude in addition to the orientation does not alter 
the qualitative behavior, we use the same function pro- 
posed in Ref. [l6| to study characteristic length scales of 
glasses, 

B(b) = (J2Uj(b) 2 ) 1/2 /d. (7) 

3 

In Eq. 0, Vj(b) denotes the coarse-grained NAD vec- 
tor Uj(6) = N^ 1 J^k dfcWb(r^), which is the analogue to 
Eq. © for the NAD instead of their unit vectors. By 
construction, the function B(b) has very similar limiting 
behavior as D{b). For large b, B vanishes since J2k=i ^fc 
corresponds to a rigid translation of the whole system 
and is conventionally set to zero. For small b, XJj — > dj 
and B approaches one. Figure |H] shows the function B(b) 
for the BMLJ model subject to homogeneous deforma- 
tion with 7 = 10 -4 . We observe a very similar behavior 
of B(b) compared to D(b). Thus, the qualitative con- 
clusions are robust and hold for different definitions of 
coarse-graining functions. For a more quantitative com- 
parison, we extract the correlation length £' B analogous 
to £' D above from the decay B(b) w exp \—b/£' R ] (see also 
the supplementary information of Ref. [30( • The growing 
length scale, extracted by least-square fitting an expo- 
nentially decaying function to the data is plotted in the 
inset of Fig. [5J The length scales £' D and £' B not only 
show a very similar temperature-dependence but are also 
quantitatively quite similar. 

A. Static length scale from histogram of 
coarse-grained NAD 

In addition to the quantities D(b) and B(b) that are 
averages of the coarse-grained NAD, we have also studied 
the distribution of the coarse-grained NAD orientations 
hb(D 2 ) and extracted a correlation length £g directly 
from them. This correlation length displays the same 
type of temperature dependence of £' B and £' D and it is 
also quantitatively consistent with them. 

For small coarse-graining distances 6, only a single 
particle contributes to the sum in Eq. (|6|). Hence the 
histogram of D 2 values is peaked around 1, hb(D 2 ) — s- 
<5(D 2 — 1). When coarse-graining over larger and larger 
distances b, Dj successively decreases and therefore hb 
accumulates more weight at small values of D 2 until hb 
becomes strongly peaked near D 2 = when b becomes 
large [3(| ■ The observed behavior bears some similarities 
to distribution of order parameters near phase transi- 
tions, where the histogram switches between mostly or- 
dered to mostly disordered states when passing through 
the transition. Here, the transition between mostly corre- 
lated (peak of hb near 1) and uncorrelated (peak near 0) 



regions happens at a value of b which strongly increases 
with decreasing temperature. Therefore, the histograms 
hb(T) 2 ) allow us to define a static correlation length £g 
as the value of b that characterizes this transition. As 
in Ref. [30| . we choose as definition the coarse-graining 
distance b for which the peak-location has shifted from 
one to 1/e. The length £b measures the average domain 
size of NADs over which displacements are correlated. 
The temperature dependence of £b for BMLJ (BMSS) 
model is shown in the inset of Fig. [TU] (Fig. ITT]) . While 
the value of £b is on the order of a particle diameter Co 
at high temperature, it increases considerably by cooling 
the system towards supercooled regions in quantitative 
agreement with the behavior of £' B and £' D discussed in 
Sec. lVII fsec also T values listed on TableUJ. Alternatively 
one also can define £b as the coarse-graining distance for 
which the variance (D 4 ) — (D 2 ) 2 is maximized. The tem- 
perature dependence of £b is robust and does not depend 
on different definitions [3(| H2 ■ 

How sensitive is the length scale on the amplitude 7 of 
the applied shear deformation? Figure IT21 shows that for 
small deformations, 7 < I0~ 4 , the correlation length £b 
becomes independent of 7. This finding allows us to in- 
terpret £b found above as an intrinsic property of the sys- 
tem, characterizing the correlation between two nearby 
inherent structure configurations. When increasing the 
deformation amplitude beyond 10~ 3 , £g decreases until 
the temperature-dependence is washed out by the shear. 
In this regime, the shear deformation is strong enough to 
decorrelate the initial and final inherent structure config- 
uration, which may be interpreted in terms of meta-basin 
transitions [53| . 

We want to stress once more that the procedure for ob- 
taining the correlation length in this work is performed 
just on simulation snapshots with static deformation. 
Hence our length scale has a purely static character and 
its significant increase upon cooling can be interpreted as 
a structural signature of different dynamical regimes in 
supercooled liquids. 



B. Finite-size scaling analysis 

Figs. [TU1 and ITT1 show that £b depends strongly on the 
linear size L of the system. The strong system-size de- 
pendence together with the significant increase of £g is 
reminiscent of critical phenomena. There, the growing 
correlation length of the order parameter characterizing 
the transition gets cut-off by the system size. To test if 
the system-size dependence of £s has any critical charac- 
ter at low temperatures, we perform a finite-size scaling 
analysis. We assume that the correlation length - as- 
sociated with an unknown order parameter for glasses - 
diverges with an exponent v at the critical temperature 
T c as £ ~ t~", where t is the reduced distance from the 
critical point, t = (T — T c )/T c . Therefore, should di- 
verge as well and in an infinite system £g ~ t~ p , where p 
is the critical exponent characterizing its divergence. On 
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the basis of the scaling hypothesis for critical phenom- 
ena [I?]], the corresponding quantity £,b,l{T) in a finite 
system of linear size L should follow the behavior 

^B.Lm^L^Q^iL^t), (8) 

where Q^ B (x) is a universal scaling function. Therefore, 
plotting £,b.l{T)L~ p / v as a function of the scaling vari- 
able tL 1 ^ should collapse all data points for different 
system sizes onto one master curve. The occurrence of 
a thermodynamic phase transition has been put forward 
by Random First Order Theory (RFOT) and related the- 
ories of the glass transition [5g, [59| and it has been dis- 
cussed whether this could take place at the temperature 
where the extrapolated configurational entropy vanishes 
Xk (the Kauzmann temperature). Fig. [TU] shows the data 
collapse obtained for the the BMLJ model using Tk as 
the critical temperature [H^, [59[ . For the BMLJ model, 
Tk is estimated numerically to be Tk ps 0.30 [36|, l37j |. 
Fixing this value for T c , we have varied the critical ex- 
ponents v and p. The best data collapse is obtained for 
the values p ps 0.9 ± 0.1 and v ps 0.65 ± 0.1, which is 
the case shown in Fig. [TU] The error bars for the critical 
exponents are estimated by varying the latter until the 
quality of data collapse starts to get worse. Because of 
very long equilibration time of the system at low tem- 
perature and of the small increase of £g for small system 
sizes, our data are still relatively far from the hypotheti- 
cal critical point and this might be one of the reasons for 
the rather large uncertainties in the values of these ex- 
ponents. The good collapse of the numerical data, how- 
ever, indicates that the critical region is large enough 
to be detected in the temperature range where we could 
equilibrate our systems. In Ref. [3(|, we have performed 
the finite-size scaling analysis by lifting the assumption 
T c = Tk- The results showed that one still could get a 
data collapse for values of T c which lay within the interval 
0.25 ^Tc^ 0.4. For the BMSS model, T K is estimated 
to be T K ps 0.11 [H|. Fixing T c to this value, the best 
data collapse is obtained by choosing the critical expo- 
nents as p ps 1.5 ± 0.2 and v ss 0.60 ± 0.15. The case is 
shown in Fig. 1111 The value of the exponent v is similar 
for BMLJ and BMSS models and this suggests that the 
divergence of the underlying correlation length £ associ- 
ated with the order parameter of the transition has the 
same origin in both cases. In Ref. [3(j, we have shown 
that also the dependence of relaxation time and configu- 
rational entropy on £b for the BMLJ model are in good 
agreement with predictions of RFOT. It is interesting to 
notice that, within the standard RFOT , [6^ , the mo- 
saic length scale £ m is expected to grow significantly only 
at very low temperatures (T < Tmct) where long-lived 
metastable states occur, whereas our results indicate that 
the characteristic length scale £,b starts to grow already 
at relatively higher temperatures. The growth of the 
so-called point-to-set correlation length at (T > Tmct) 
has also been observed in recent numerical simulations 
of three glass formers in Ref. [6l|. Overall, the picture 
emerging from our results is qualitatively consistent with 



the RFOT scenario. On the other hand, the uncertainties 
in the values of the critical exponents, as well as differ- 
ent results obtained in other recent studies |47|, indicate 
that further quantitative investigations are still needed to 
better clarify the nature of the critical behavior observed 
here. 



VII. CONCLUSIONS 

Cooperatively rearranging regions (CRR) play an im- 
portant role for the dynamics of supercooled liquids with 
their sizes growing moderately upon approaching the 
glass transition (j-Q • Very recently, structural signatures 
of the CRR have been detected, suggesting medium- 
range order Q or localized soft modes Q as triggers of 
the rearrangements. Here, we have shown by extensive 
molecular simulations that correlations in neighboring IS 
are quite reminiscent of CRR that are observed in the 
system's dynamics [||. The behavior we report here is 
found for different models of fragile glasses, i.e., the Kob- 
Andersen model and a binary soft sphere mixture. The 
distance between two IS - that are related via shear defor- 
mation of rather small amplitude 7 - sharply decreases 
with decreasing temperature below the onset tempera- 
ture of the landscape-dominated regime ■ We observe 
that the mean distance between the two IS varies expo- 
nentially with the inherent structure energy eis- Further- 
more, we observe a crossover between two regimes that 
is also present in a qualitative change of the NAD distri- 
bution from an exponential to a power-law shape. The 
exponent of the latter can be rationalized from elasticity 
arguments [49}. We use different measures for correla- 
tions of NAD vectors. Qualitatively, we find similar re- 
sults as experiments on correlated motions in colloids . 
Quantitatively, the static correlation length extracted 
from NAD of IS grows significantly upon getting closer to 
the glass transition. Moreover, the finite-size scaling of 
our results indicates that the correlation length diverges 
at low temperatures. Since thermal fluctuations tend to 
wash out the correlations in the NAD, our method is very 
sensitive and provides an efficient tool to investigate cor- 
related regions. This is a first interesting outcome of 
this work. It is worth noting, with this respect, that 
NAD field has proven to be an insightful investigation 
tool also for the mechanics of amorphous solids [16[ and 
hence the approach we propose has a good potential to 
bridge the investigation of supercooled liquids to the one 
of amorphous solids, within a unifying picture of glass 
transition. In addition, the results here discussed sug- 
gest that the long range spatial correlations detected by 
the NAD field upon lowering the temperature might be 
indeed the spatial correlations underlying the onset and 
development of cooperative dynamics typically observed 
in supercooled liquids. This would be a significant new 
insight into the physics of supercooled liquids and more 
work is currently in progress to clarify this issue. The 
critical growth of spatial correlations that we have re- 
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ported here overall supports the glass transition scenario 
based on Random First Order Theory and the analy- 
sis done seems to be quite robust. On the other hand, 
the finite-size scaling proposed also raises the question of 
identifying the critical temperature and for the moment 
we cannot rule out different possibilities: this requires 
further and more quantitative analysis. Finally, it would 
be extremely interesting, at this point, to be able to di- 
rectly relate the growth of the correlation length and the 
features of the NAD field to the change in the transport 
properties of the system, i.e. in its viscoelastic response: 
whereas, to some extent, we can relate our results to the 
viscosity increase through the RFOT [30| , a more coher- 



ent connection with the mechanical response and with 
the onset of rigidity in the material needs to be devel- 
oped. 
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FIG. 2. Main frame: Average bond orientational order mea- 
sured in the starting configuration X and its inherent sruc- 
ture X q as a function of temperature for BMLJ model. Inset: 
Average energy of the inherent structure configuration as a 
function of temperature T for BMLJ model. Error bars are 
smaller than the size of the symbols. 
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FIG. 3. Distribution of mismatch components hi (defined 
in the text) in the BMLJ model for homogeneous deforma- 
tion of strength 7 = 10~ 4 . Temperature is increasing as 
T = 0.43,0.55,0.66,1.0 from left to right at the bottom of 
the figure. At high temperatures the distributions are expo- 
nential (dashed line shows a fit exp [— |d a |/0.1] to the data at 
T = 1.0), but start to develop a power-law tail at low tem- 
peratures. The solid line shows a fit to a power-law with an 
exponent —5/2. 




FIG. 4. NAD field for for one selected configuration at high 
(T = 1.0 left) and low (T = 0.42 right) temperature. Ho- 
mogeneous deformation of strength 7 = 10~ 4 was applied to 
BMLJ system. For clarity just particles on visible surfaces 
are shown. 
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FIG. 5. The mean length of the mismatch vectors d defined 
in the text as a function of temperature for the BMLJ model 
and for different deformation magnitude 7. 
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FIG. 6. Main frame: The mean length of NAD vectors as a 
function of inherent structure energy eis for the BMLJ model. 
Note the crossover from high temperature to low temperature 
regime which happens slightly below T a . Inset: temperature 
dependence of d and eis. 
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FIG. 7. Correlation of longitudinal mismatch directions C\\ 
as a function of particle separation r defined in Eq. ([5]). Ho- 
mogeneous shear deformation with 7 = 1CP 4 was applied to 
BMLJ model with N = 8000 particles. Inset shows the length 
scale obtained by fitting an exponentially decaying function 
to data. 
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FIG. 8. (D 2 ) 1/2 defined in Eq. © as function of b for selected 
temperatures T. Results are shown for the BMLJ model un- 
der homogeneous shear with 7 = 10 -4 . Inset shows the grow- 
ing length scale obtained by fitting an exponential decaying 
function to data points which are shown in main figure. 
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FIG. 9. The coarse-graining function B(b) denned in Eq. Q 
for selected temperatures T. Results are shown for the BMLJ 
model under homogeneous shear with 7 = 10 -4 . Inset shows 
the growing length scale obtained by fitting an exponential 
decaying function to data points which are shown in main 
figure. 




FIG. 10. Data collapse for the BMLJ model. Inset: Static cor- 
relation length £b, extracted from histogram of coarse-grained 
NAD orientations as function of temperature for different sys- 
tem sizes. 
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FIG. 12. The length scale extracted from peak-location of 
histograms of coarse-grained NAD orientations as function 
of temperature T for different strength of shear deformation 
7 for BMSS model with N = 16000. At low 7 the length 
becomes independent of deformation magnitude and hence 
an intrinsic property of the system. 



